Cluster update for tensor network states 
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We propose a novel recursive way of updating the tensors in projected entangled pair states by 
evolving the tensor in imaginary time evolution on clusters of different sizes. This generalizes the so- 
called simple update method of Jiang et al. [Phys. Rev. Lett. 101, 090603 (2008)] and the updating 
schemes in the single layer picture of Pizorn et al. [Phys. Rev. A 83, 052321 (2011)]. A finite-size 
scahng of the observables as a function of the cluster size provides a remarkable improvement in the 
accuracy as compared to the simple update scheme. We benchmark our results on the hand of the 
spin 1/2 staggered dimerized antiferromagnetic model on the square lattice, and accurate results for 
the magnetization and the critical exponents are determined. 
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Numerical simulation of strongly correlated quantum 
systems in dimension above 1 remains one of the big chal- 
lenges in condensed matter physics. Quantum monte 
carlo fails for models suffering from the sign problem, 
and the density matrix renormalization group method 
(DMRG) has a problem in 2D because the violation of 
the area law of entanglement entropy [l[. Tensor net- 
work states (TNSs) 0-0] ' the other hand, naturally 
generalize MPSs to higher dimensions and fulfill the area 
law of entanglement entropy, and are proving to provide 
powerful tools to simulate problems in above 1 dimen- 
sion [5-16|. 

One of the major difficulties in a TNS algorithm is the 
scaling of the computational demands as a function of 
the virtual bond dimension D. Inspired by the DMRG 
algorithm, where truncation is made with respect to the 
reduced density matrix, a contraction al gori thm at the 
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wavefunction level was recently proposed 
the rapidly growing renormalized physical index with the 
system size causes a barrier that hampers the efficiency 
of the algorithm. The simple update [181], proposed as a 
generalization of the infinite time-evolving block decima- 
tion (iTEBD) algorithm ^ to TNSs, successfully avoids 
the exponentially large Hilbert space, and is very effi- 
cient [i3,[ii,[2£ 



2l| . However, the product environment 



is too simple to capture the long range entanglement and 
correlations near the critical point of a second order phase 
transition 13j, |2l| . The efficiency of the simple update 



still sheds light at controlling the Hilbert space in the 
complete contraction algorithm jlTj . The goal of this pa- 
per is to demonstrate that those methods can be merged 
together in such a way that the advantages of both meth- 
ods are preserved. 

Method - The wavefunction of an infinite projected en- 
tangled pair state (iPEPS) in terms of local tensors on 2 
sites in the simple update [l^ is given by 
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where \se) = \h jTk,m,n,p) represents the environmen- 
tal degrees of freedom. We generalize this construction 
to a block of x ly sites with the open virtual bonds 
weighted by the diagonal tensors \/As, where As are 
the entanglement spectra with respect to each bond, 
as in Fig. [TJa). As a consequence, the infinite Hilbert 
space is reduced to the product space of the open virtual 
bonds and the remaining physical bonds in the cluster as 
Fig.[TJb). The long range entanglement is gradually con- 
sidered by taking larger block sizes. Due to the modifica- 
tion of the boundary sites of the cluster, an iPEPS wave- 
function is converted to an open boundary (OB) PEPS 
of size Ix X ly, of which the contraction can be done effi- 
ciently at the wavefunction level. This conversion is only 
made at the evolution stage. Once the ground state ten- 
sors are obtained, one evaluates expectation values such 
as the energy and rnagnetization using the Monte Carlo 
sampling technique |2C| to a system with periodic bound- 
ary (PB) condition and of system size L x L. 
The iterative steps of the cluster update are: 

• write the iPEPS wavefunction into an OB PEPS us- 
ing the As from the previous iterations as Fig.[IJa), 

• act with a local (imaginary time) evolution oper- 
ator on several sites located in the center of the 
Ix X ly cluster as Fig.[TJb), 

• contract every virtual bond except the one with 
the evolution operator, obtain the projectors and a 
new A as discussed below, and project the enlarged 
bond back to D, 

• finally replace the iPEPS wavefunction by the up- 
dated tensors as Fig. [Ijc), and move to a different 
bond. 

There are two main choices to be made in the above 
update procedure: first of all, how to contract the en- 
vironment tensor efficiently to two sites under the time 
evolution; second, how to obtain the updated tensor. 

To answer the first question, let us remind the conven- 
tional way of contracting the environment tensor. In the 
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FIG. 1: (a) Wavefunction of an infinite lattice in terms of 
Ix X ly local tensors and the entanglement spectra placed on 
each open bond, (b) An equivalent OB PEPS to (a), with 2 
sites acted by evolution operator in the center of the cluster, 
(c) After projection with projectors Pa, Ps, and replace A, B 
with A, B everywhere, a new iPEPS wavefunction is obtained. 
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FIG. 2: (a) MPSs \pb) and \pt) of a conventional contraction 
method dealing with the norm of a wavefunction 

(b) MPSs \(j)t) and l^j) of a novel contraction method dealing 
with the wavefunction lip) directly. 
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FIG. 3: (a) Renormalize a row into the MPS {(pb). (b) Ab- 
sorb a row into \(j)f,) enlarges both the internal bond and the 
physical bond by a factor D and d respectively, (c) Projectors 
(blue triangles) PJ^, P^ and P5 is calculated to reduce both 
the internal bond and physical bond of matrix M'. 
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FIG. 4: (a-d) Calculate the residue matrix (blue diamonds) 
Ri (Li) of the left (right) half chain up to site i by iteratively 
using QR (LQ) decompositions, (e) Replace left (right) half 
chain by (L'^^) and use an SVD to reduce the dimension 
of the j-th bond, (f) Inserted an identity (hollow diamond) 



at bond i in original MPS and replace R'L 
obtain the projectors P^ and P^^. 
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conventional evolution methods, one constructs a double 
tensor network, which represents the norm of the 

state lip), then contracts this norm to obtain the envi- 
ronment tensor of two sites. The major approximation 
is to represent the bottom (top) rows of the norm by 
an MPS: \pbilj)) = TTk{Ml~k)\MjM i\pt(i'J)) = 
Tr^,{Mi',k')\Mf,k'))) as Fig. [^a), where \Mh~k)) 
(|0f (i', fc'))) is the wavefunction of the bottom (top) 
rows with the highest (lowest) vertical bonds open, as 
Fig. EJb), i = {ii, 12, *3, *4}, etc. However, during the 
contraction, one ignores the fact that \ph(i^j)) is hermi- 
tian if written as a density matrix of the fictitious degrees 
of freedom i and j, and treats \pb{i,j)) as a general MPS 
with degrees of freedom {ii,ji} at each site /. This will 
cause ill-condition if the truncation error is large. We 
therefore contract the wavefunction directly to ensure the 
hermiticity of the norm. The degrees of freedom at each 
site I of \4>b{i, k)) is {i;, fc;} with ii represents the topmost 
vertical index and fc; represents the renormalized physi- 
cal index of the bottom-half-column. The argument 
for {(pbih k)) equally applies to |0t(j', k')) throughout the 
paper. 

To contract a tensor network wavefunction with OB 
conditions, one absorbs a row into \4>b) to obtain \(j)b), as 
Fig. [2Ja,b). Unlike in the case of contracting the norm, 
this is described as a matrix product operator (MPO) 



acting on an MPS, in such a way that not only the in- 
ternal bond dimension x but also the physical degrees 
of freedom d of \(j>b) is increased by a factor of D and d 
respectively as Fig. [3{b), where D and d are the virtual 
bond dimension and the physical degrees of freedom in 
the original lattice. In order to keep the renormalization 
under control, one has to truncate both, as Fig.[3l^c). 

To reduce the bond dimension between site i and i + 1 
of an OB MPS, a singular value decomposition (SVD) is 
done on \(j)b). As in Fig. [Sje), one needs the right (left) 
residue matrix R* (L*"*"^) of the left (right) half chain 
up to the site i (i -I- 1) to bring them into their normal 
form [ij. To be specific, one calculates the boundary 
residue matrix (L") using a QR (LQ) decomposition, 
n is the length of \(t>b), 

r' 

then moves one site to the right (left) as in Fig. IDJa-d), 

I' r' 
r' /' 
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FIG. 5: (a) Reduced density matrix of the physical bond s at 
site i of \(t>b)- (b,c) A reformation of (a) by M' and the left 
and right residue matrices L*"*"^, R'~^. 



where is the matrix at the site i of \4>b), l,r,u, s rep- 
resent the_left, right, up and physical index respectively, 
and Q' (Q*) is isometric matrix. Thus the left (right) 
half chain, which is replaced by (L'+^), is in its nor- 
mal form as^Fig. S^e). An SVD R^L^+i w IJSVt will 
minimize \\(l>b) — |0f,)|Pi where is the one-bond re- 
duced MPS to \<j>b), bar means taking the leading sin- 
gular values or vectors. To derive the projectors for this 
reduction, one inserts an identity (R*)~-'^R*L*+-'^ (L*+-'^)^-'^ 
at bond i in the original MPS and substitutes R*L'+^ by 
USVt as Fig. Iljf). Equally distributing S to each side, 
one writes the projectors at bond i as 



To avoid the matrix inversion, one can use (R*L*+^)^^ 
V^Ut to rewrite the projectors as 



(6) 
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note that taking will not cause a singularity as the 
small singular values are discarded. To reduce the phys- 
ical degrees of freedom of M', one considers the j;educed 
density matrix pi ,,, of site i as Fig. [5l Define M^^^^ — 

Eur'Kl'^i'usr'KXl, then Pi,, = ElurMliuMl^s'- 

The leading eigenvectors Pg of p* = PgA*Pg is the pro- 
jector to reduce the physical degrees of freedom of M*. 
Upon obtaining P|^, P]^ and Pg, one projects according 
to Fig. |3l[c). Note that one left (right) sweep can obtain 
R* (L*) for all sites. It turns out that the reduction of 
each internal bond and physical bond can be done simul- 
taneously as if they are independent. This concludes the 
question of the contraction at the single-layer wavefunc- 
tion level. 

Now we focus on the question of how to obtain the 
evolved tensors using the partially contracted wavefunc- 
tion. After renormalizing rows from below (above) into 
{(f'b) (I </**)) as Fig. EJa), one renormalizes columns from 
left (right) into \(f)i(h,si)) {\(j)r(h,s ,■))), as in Fig. [Hb), 
where the thickened bond is due to the action of the evo- 
lution operator. Define 



\tp{si,Sr)} = ^ |0i(/l. Si)) ® (f)r{h, Sr)), 



(8) 



where h = {hi, h2, h^} is the horizontal internal bond, 
the projector to reduce the thickened bond has to mini- 
mize the 2-norm \ \i/j{si , Sr)) — \ {si , Sr))]]"^ ■ Unlike in the 
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FIG. 6: (a) Bottom and top MPS \<j)b), \(t>t) sandwich a row 
with the evolution operator, (b) Further absorb all columns 
into the left (right) MPS (|0r)). (c) Gontract all virtual 
bonds in (b) except h2 to form a single tensor T. 



OB MPS, where one cut will bipartite the wavefunction, 
here one deals with a ring MPS. We present an empirical 
way to calculate the projectors. First contract all virtual 
bonds in {ipisi, Sr)) except the thickened bond /12 to form 
a single tensor T as Fig.[6l^c); in analogy to an OB MPS, 
make a QR (LQ) decomposition to calculate the right 
(left) residue matrix R (L) of the tensor T as 



Trs,l = Qrs.l'Rl'.l, 
I' 

Tr^sl = ^ ' ^r.r'Qr' .sh 



(9) 



(10) 



where Q (Q) is isometric matrix; insert R^^RLL"^ be- 
tween the l,r indices of tensor T, in analogy to Eq. [71 
derive the projectors to the imaginary time evolution as 

Pa=LV^, (Pb)^ = 4^U^R, (11) 
V A V A 

where an SVD RL « UAVt is performed. A for this 
particular bond contains the entanglement spectrum that 
replaces the previous one for the next evolution step, 
analogous to what happens in the simple update scheme. 
Calculating the T tensor explicitly is expensive, what 
one does instead when dealing with the object Fig. [6jc) 
is treating the r index as the physical ones {si,Sr} to 
calculate the R matrix and again treating the I index as 
the physical ones to calculate the L matrix. Another way 
to update the projectors Pa, Pb would be alternatively 
solving the multi-quadratic equations [l[, where the T 
tensor is explicitly needed. 

Results - We benchmark our method using the spin 1/2 
staggered dimerized antiferromagnetic Heisenberg model 
on square lattice with PB conditions and compare our re- 
sults to that obtained from the simple update ^21,] , iPEPS 
method ^ and the QMC simulation The Hamil- 
tonian is written as 



H 



(id) 



J' ^ 

{(fe,/» 



S,, (J,J'>0), (12) 



where ((fc, 0) a-re the horizontal nearest neighbor pairs 
satisfying mod(a;fe,2) = mod(yfe,2), Xk,yk is the lattice 
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FIG. 7: (a) Ground state energy errors per site for L — 16 us- 
ing tensors of bond dimension D = 5 obtained via the simple 
update {Ix = 2) and the cluster update {Ix ~ 4, 6, 8). 
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FIG. 8: (a) Sublattice magnetization square for L = 16 using 
the same parameters as Fig. [T] black solid line is the SSE 
results for L = 16, solid circles are the CSE results from (b). 
(b) The cluster size extrapolations (CSE) for the finite size 
{L = 16) magnetizations using cluster sizes Ix ~ 2, 4, 6, 8. 



coordinate of site fc, {i,j) are all other nearest neighbor 
pairs, and J, J' are the coupling strengths. This model 
is frustration free and has been extensively studied 12311 
using the stochastic series expansion (SSE) method [2J]. 
Increasing the coupling strength ratio J'/J drives the 
system through a second order phase transition from the 
antiferromagnetic ordered phase to a magnetically disor- 
dered phase. We set J = 1 for convenience hereafter. 

We obtain the finite size ground state energies and the 
magnetization for systems with linear size L = 4, 8, 16. 
Each system is measured using the variational monte 




FIG. 9: (a) Finite size magnetizations and their thermody- 
namic limits via the ESS from (b), the black solid lines are 
the SSE simulation results for system size L = 4, 8, 16 and the 
black dash line is their FFS results, (b) Finite size scaling of 
the magnetization square using system sizes L = 4, 8, 16. 
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FIG. 10: Fitting of the magnetization curve for the critical 
value Jc and the critical exponent /3. 



carlo (VMC) sampling technique [2(| taking the tensors 
of bond dimension D = 5 obtained through the cluster 
update with cluster SIZG X ly — 2 X 1, 4 X 4, 6 X 6 
and 8x8. For this specific model, the simple update is 
equivalent to a cluster update with cluster size 2x1. 

We present the ground state energy errors per site for 
system size L = 16 with D = 5 m Fig. [7] and the sub- 
lattice magnetization square defined as following f23| in 
Fig.Ii 
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(13) 



Significant improvement in energy and magnetization is 
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achieved by increasing the cluster size Ix from 2 to 8. We 
make a cluster size extrapolation (CSE) for the finite size 
magnetization in Fig. |HJb) as 



(14) 



M^J',L,Q = M^J',L) + aj,/^l, 



where M^{J',L) and aj' are the fitting parameters. The 
CSE results are in good agreement with the finite size 
magnetization through the SSE simulation, as plotted in 
Fig. IHl^a). Here the cluster size plays the role of the 
bond dimension D to provide a scaling scheme. To obtain 
the sublattice magnetization in the thermodynamic limit, 
we use the following finite size scaling (ESS) formula [23- 
^ to extrapolate M (J') as in Fig. [^b) 



M^{J',L) = M^{J') + bj,/L, 



(15) 



where M'^{J') and bji are the fitting parameters. This 
FSS formula is only valid for the states with persisting 
sublattice magnetization, thus the extrapolated negative 
M^{J') simply means that these states have no magnetic 
order. The finite size magnetizations M{J',L) together 
with their thermodynamic limits M{J') are plotted in 
Fig.lHIa). The extrapolated Af(J') lies a bit off the SSE 
FSS results near the critical point, because the latter are 
fitted with a sub- leading corrections 1/L^ using system 
sizes L = 4, 8, 16, 32. To determine the critical value 
and the critical exponent /3, we fit M{J') by 



M( J') = c( ~ J') 



/3 



(16) 



where J^, /3 and c are the fitting parameters, as in 
Fig. [TUl We obtain a critical value = 2.501(1) and 
the critical exponent (5 = 0.37(1). The SSE results are 
= 2.5198(3) and /3 = 0.376(5) 23]. The offset in 
is due to ignoring the sub-leading corrections in the 
FSS formula Eq. 1151 Considering it would require results 
from a much larger bond dimension D for system size 
L = 32. The critical value from the simple update is 
= 2.56 and from the iPEPS is « 2.85 [H. 
discussion - The cluster update allows one to accu- 
rately determine the behavior of a second order phase 
transition using a tensor network state of intermediate 
bond dimension D and relatively small cluster size l^. 
The simple update is efficient, however it always gener- 
ates a fat tail near the exact transition point; on the con- 
trary, the cluster update that has been introduced here 
improves those results significantly. This improvement is 
especially important when a narrow intermediate phase 
is present in the phase diagram, such as in the frustrated 
Ji — J2 Heisenberg model on square lattice. Comparing to 
the complete contraction algorithm in 17[, the Hilbert 
space is drastically reduced, which significantly boosts 
the efficiency, and the algorithm is simpler in the sense 
that it is free of dealing with the boundary evolutions. 
The cluster update scales as D^x^dd. The renormalized 
physical index d varies depending on the entanglement 



of the state, however it is bounded by D~x^ . Choosing a 
small cluster size also relaxes the scaling of d, e.g. the 
complexity with cluster size 2x2 scales as D^, and 4x4 
scales as if taking x = Z?. 

Conclusion - We presented a cluster imaginary time 
evolution method for a tensor network state (TNS) de- 
scribing the ground state of strongly correlated quantum 
systems. We benchmarked this method with the stag- 
gered dimerized antiferromagnetic Heisenberg model on 
the square lattice and accurately determined its critical 
value and critical exponents (3 using a TNS with fairly 
small bond dimension Z? = 5; this provides clear evidence 
for an improvement over the simple update scheme in 
[3] ■ The efficiency and accuracy of this method should 
allow tensor network simulations to be applied to a zoo of 
interesting models that are not easily accessible by other 
methods, especially as large values of D can be treated 
[23. 

Acknowledgments: This project is supported by the 
EU Strep project QUEVADIS, the ERC grant QUERG, 
and the FWF SFB grants FoQuS and ViCoM. The com- 
putational results presented have been achieved using the 
Vienna Scientific Cluster (VSC). 



[1] F. Verstraete, V. Murg, and J. I. Cirac, Adv. Phys. 57, 
143 (2008). 

[2] N. Maeshima, Y. Hieida, Y. Akutsu, T. Nishino, and 
K. Okunishi, Phys. Rev. E 64, 016705 (2001). 

[3] F. Verstraete and J. I. Cirac, arXiv:cond-mat/0407066vl. 

[4] G. Vidal, Phys. Rev. Lett. 101, 110501 (2008). 

[5] V. Murg, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 
95, 057206 (2005); V. Murg, F. Verstraete, and 
J. I. Cirac, Phys. Rev. A 75, 033605 (2007); V. Murg, 

F. Verstraete, and J. I. Cirac, Phys. Rev. B 79, 195119 
(2009); V. Murg, 6. Legeza, R. M. Noack, and F. Ver- 
straete, Phys. Rev. B 82, 205105 (2010). 

[6] R. Orus, A. C. Doherty, and G. Vidal, Phys. Rev. Lett. 
102, 077203 (2009); R. Orus and G. Vidal, Phys. Rev. B. 
80, 094403 (2009). 

[7] J. Jordan, R. Orus, G. Vidal, F. Verstraete, and 
J. I. Cirac, Phys. Rev. Lett. 101, 250602 (2008); J. Jor- 
dan, R. Orus, and G. Vidal, Phys. Rev. B 79, 174515 

(2009) . 

[8] G. Evenbly and G. Vidal, Phys. Rev. Lett. 104, 187203 

(2010) . 

[9] P. Corboz, J. Jordan, and G. Vidal, Phys. Rev. B 82, 
245119 (2010); P. Corboz, S. R. White, G. Vidal, and 
M. Troyer, Phys. Rev. B 84, 041108 (2011); P. Cor- 
boz, A. M. LiiucMi, K. Penc, M. Troyer, and F. Mila, 
,arXiv:11 08.2857 
[10] Z.-C. Gu, M. Levin, and X.-G. Wen, Phys. Rev. B 78, 
205116 (2008); Z.-C. Gu, arXiv:1109,4470; Z.-C. Gu, 
H.-C. Jiang, D. N. Sheng, H. Yao, L. Balents, and X.- 

G. Wen, arXiv:1110.1183 

[11] Q.-Q. Shi, S.-H. Li, J.-H. Zhao, and H.-Q. Zhou, 
arXiv:0907.5520; S.-H. Li, Q.-Q. Shi, and H.-Q. Zhou, 
,arXiv:1001.3343. 



6 



[12] H.H. Zhao, Z.Y. Xie, Q.N. Chen, Z.C. Wei, J.W. Cai, and 
T. Xiang, Phys. Rev. B. 81, 174411 (2010); H. H. Zhao, 
C. Xu, Q. N. Chen, Z. C. Wei, M. P. Qin, G. M. Zhang, 
and X. Tao, arXiv: 1105.2716 

[13] W. Li, S. S. Gong, Y. Zhao, and G. Su, Phys. Rev. B 81, 
184427 (2010); W. Li, S. S. Gong, Y. Zhao, S. J. Ran, 
S. Gao, and G. Su, Phys. Rev. B 82, 214431 (2010). 

[14] L. Wang, Y.-J. Kao, and A. W. Sandvik, Phys. Rev. E 
83, 056703 ( 2011); J.-F. Yu and Y.-J. Kao, 
larXiv:1108.3393l 

[15] C. Liu, L. Wang, A. W. Sandvik, Y.-C. Su, and Y.- 
J. Kao, Phys. Rev. B 82, 060410R (2010); C. Liu, D.- 
X. Yao, and A. W. Sandvik, arXiv:1110.0761 

[16] C. V. Kraus, N. Schuch, F. Verstraete, and J. I. Cirac, 
Phys. Rev. A 81, 052338 (2010); P. Corboz, 
R. Orus, B. Bauer, and G. Vidal, Phys. Rev. B 
81, 165104 (2010); Z.-C. Gu, F. Verstraete, and X.- 
G. Wen, arXiv:1004.2563, I. Pizorn and F. Verstraete, 
Phys. Rev. B 81, 245110 (2010). 



[17] I. Pizorn, L. Wang, and F. Verstraete, Phys. Rev. A 83, 
052321 (2011). 

[18] H. C. Jiang, Z. Y. Weng, and T. Xiang, Phys. Rev. Lett. 

101, 090603 (2008). 
[19] G. Vidal, Phys. Rev. Lett. 98, 070201 (2007). 
[20] L. Wang, L Pizorn, and F. Verstraete, Phys. Rev. B 83, 

134421 (2011). 

[21] C.-Y. Huang and F.-L. Lin, Phys. Rev. B 84, 125110 
(2011). 

[22] B. Bauer, G. Vidal, and M. Troyer, J. Stat. Mech. 2009, 
P09006. 

[23] S. Wenzel, L. Bogacz, and W. Janke, Phys. Rev. Lett. 

101, 127202 (2008). 
[24] A. W. Sandvik, Phys. Rev.B 56, 11678 (1997). 
[25] H. Neuberger and T. Ziman, Phys. Rev. B 39, 2608 

(1989). 

[26] D. S. Fisher, Phys. Rev. B 39, 11783 (1989). 
[27] L. Wang et al.. Work in Progress. 



